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OBJECTIVE 


Improve  on  existing  ray-tracing  methods  used  in  minicomputer-based  propagation 
assessment  systems.  Two  problems  are  addressed:  (I)  improvement  in  accuracy  and  (2) 
speed  of  calculations. 


RESULTS 

Exact  ray  paths  can  be  calculated  rapidly  with  the  equations  in  this  report.  They 
are  based  on  two  assumptions  which  should  be  considered  in  their  use,  however.  These  are 
that  (1)  the  earth’s  magnetic  field  and  (2)  any  spatial  variation  in  the  ionosphere  profile  can 
be  ignored. 


RECOMMENDATIONS 

The  ray-tracing  equations  presented  here  should  be  considered  for  implementation 
in  mini-  and  microcomputer-based  hf  assessment  systems.  In  such  implementation,  the 
ionosphere  profile  should  be  changed  for  each  ray  hop  according  to  known  variations  in  the 
ionosphere  along  the  hf  circuit. 
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INTRODUCTION 


I 


In  hf  communication  forecasting,  it  is  often  convenient  to  use  simplified  ionospheric 
models  to  determine  the  geographical  coverage  of  radio  waves.  The  process  of  determining 
the  radio-wave  path  in  the  ionosphere  is  called  ray  tracing.  Ray  tracing  requires  computer 
calculations  which  arc  very  extensive  if  the  exact  details  of  the  electron  density  distribution 
and  the  earth’s  magnetic  field  are  included.  For  some  applications  it  suffices  to  ignore  elec- 
tron collisions,  magnetic  field  effects,  and  precise  electron  density  profile  details.  It  will  be 
convenient  to  assume  the  vertical  profile  to  be  homogeneous  along  the  path  between  the 
transmitter  and  receiver. 

The  parabolic  ionosphere  is  a popular  approximation  to  the  electron  density  distribu- 
tion in  both  the  E region  and  the  F region.  Croft  and  Hoogasian  (reference  1)  have  shown 
that  ray-path  integrals  can  be  evaluated  exactly  in  closed  form  if  the  electron  distribution  is 
modified  slightly  from  a parabola  to  what  is  called  a quasiparabola: 
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where  Ne  = electron  density;  Nm  = maximum  value  of  Ne;  r = radial  distance  from  earth’s 
center;  rm  = value  of  r where  Ne  = Nm;  r^  = value  of  r at  layer  base;  and  ym  = rm  - rb,  the 
layer  semithickness.  The  quasiparabola  differs  from  the  parabola  by  the  multiplier  (rb/r)-. 
This  factor  is  very  nearly  unity  in  the  layer  so  that,  for  practical  purposes,  the  quasiparabola 
is  indistinguishable  from  a parabola.  Its  advantage  arises  in  the  solution  of  the  ray  equations. 

Reference  1 gives  equations  for  three  ray-path  variables:  D,  the  distance  traversed, 
measured  along  the  earth’s  surface;  P',  the  group-path  distance  (signal  transmit  time  multi- 
plied by  c);  and  P,  the  phase-path  distance  (the  wave-front  transmit  time  multiplied  by  c). 

In  this  report,  additional  equations  are  presented  for  the  ray-path  coordinates  along 
the  path.  General  ray  paths  are  considered,  including  transit  through  the  layer  or  through 
partial  layer  segments.  Complicated  multisegmented  layer  profile  calculations  are  illustrated. 
Ray  paths  trapped  in  a valley  layer  are  determined  along  with  “whispering  gallery”  conditions. 


RAY-PATH  EQUATIONS 

» Reference  1 shows  that  the  integrals  for  D can  be  evaluated  in  the  following  manner 
(see  figure  1): 
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: dr 
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1.  Croft,  TA,  and  H Hoogasian,  Exact  Ray  Calculations  in  a Quasiparabolic  Ionosphere,  Radio  Science  3 
(New  Series),  No  1 , p 69-74,  1968. 
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Figure  1.  Ray-path  geometry. 

where  the  refractive  index,  n , is  defined  by 

, 80.62  Ne  - , 

V1  = 1 - — 2 = 1 - (fc/02  + [(rm  - r)fcrb/(ymfr)l  2 (3) 

and 

13  = angle  of  ray  path  to  the  horizontal 
p0  = Patr=rQ 
f = operating  frequency 

fc  = critical  frequency  of  the  layer  (f2  = 80.62  Nm)  MKS  units 
Tq  = earth  radius  (6371  km) 
rt  = r at  top  of  ray  path. 

The  ray  path  is  a straight  line  in  the  region  Tq  < r < rj,,  which  is  free  space  with  n = 1 . 
In  the  ionospheric  portion,  r > rb  (ie,  substitution  of  equation  3 into  equation  2),  the  radical 
becomes 


(4) 


r2//2 *  - 1q  cos2  Pq  = Ar2  + Br  + C 

with 

a = i - (fc/o2 + (fcrb/f  ym)2 

B 2rm^crb^  ^01^" 

c ■ (fcWf  ym)2  - ro cos20  o- 

The  coordinates  (Dr,  r)  of  a point  on  the  ray  when  r < rj,  are  on  one  of  two  straight 
[ rQ(/3  - 0Q)  (upgoing  line) 

(downgoing  line) 

where  cos  0 = (rg/r)  cos  0q.  When  the  ray  is  in  the  layer. 


lines 


Dr  = 


(5) 


D-r0<P-fio) 


^ „ rmrb 

ru  < r < 

b rb-Vm 

and  equation  2 becomes 


r 

D,-'0<t-t0>  + r0a*t0  f ~~  • 

JT  r vAr“  + Br  + C 


0 

The  integral  can  be  evaluated  by  means  of  standard  forms  given  in  many  tables  (eg,  refer- 
ence 2).  The  result  is 


ro cos  0O  \ r<2c + rbB + l 

»r"6<ei.-w+-^-  8"{rb(2C;rB;2ygr7  f 


iv/exT  J 


(6) 


where  X = Ar2  + Br  + C,  X^,  = r^  - Tq  cos2  0q,  and  0b  = cos  * ( Tq/t b cos  0q). 

Three  cases  must  be  considered,  depending  on  the  values  of  f and  0.  If  f and  0 are 
large  enough,  the  ray  will  not  return  from  the  ionosphere  to  the  ground  but  will  proceed 
through  to  the  free  space  above  the  layer.  The  direction  taken  by  the  ray  path  is  indicated 
by  the  roots  of  the  quadratic  equation  Ar2  + Br  + C = 0.  If  the  roots  are  complex 


2.  Hill,  JR,  An  Improved  Algorithm  Relating  the  F-Layer  Peak  to  M(3000)F2,  Union  Radiologique 

Scientifique  Internationale,  Boulder,  CO,  1975.  Text  available  in  NELC  TN  3097,  Naval  Ocean 

Systems  Center,  San  Diego,  CA  92152. 
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(B“  < 4AC),  the  ray  penetrates  the  layer.  When  B“  > 4AC,  the  ray  returns  to  the  ground 
after  reaching  a maximum  height,  hj  = rj  - Tq,  where 


B + ^/B  - 4AC 
2A 


When  r = rt,  equation  6 simplifies  to 


Dt  = rO(0b-0O)  + 


Tq  COS  Pq  f(2C  + Brb  + 2x/o^) 


rb\/B*  - 4AC 


It  should  be  noted  that  D(  is  just  half  the  total  ray-path  distance,  since  it  corresponds  to  the 
upgoing  portion.  Thus,  D = 2Df. 

The  third  case  is  B-  = 4AC,  which  produces  the  so-called  Pedersen  ray,  for  which  r 
approaches  r^  asymptotically.  The  launch  angle,  j3p,  of  the  Pedersen  ray  is  useful  since  it 
determines  the  angle  of  a vertical  cone  centered  at  the  antenna  above  which  rays  penetrate 
the  ionosphere.  An  optimum  antenna  will  beam  energy  below  the  cone  to  achieve  best 
surface-to-surface  communication.  The  Pedersen  ray  angle  is  found  by  calculating  the 
maximum  in  r^ 


rt,max  = "B/2A  • 

Next,  solve  for0p  from  Ar“  max  + Brj  max  + C = 0 and  C from  equation  4.  The  result  is 


rn  cos  Pn  =\/-B(rm  + B/2A1/2  . 


Dr  changes  rapidly  as  a function  of  r when  the  ray  is  near  rb  It  is  convenient  to 
invert  equation  6 to  obtain  an  equation  for  r as  a function  of  Dr.  The  result  is 


where 


(V  - B)“  - 4AC 


2C  + rbB  + 2y/cX^ 


y/c{Dr-r0(Ph-P0)l 

2 „ 
rQ  cos  PQ 


The  limits  on  Df  are  such  that 


rK  <r< 


rb-yr 


and  can  be  determined  by  equation  6. 


GROUP-PATH  DISTANCE 


The  equation  for  the  group-path  distance  can  be  derived  in  a similar  manner. 


V = r ds  = C r dr 

T~J  » " / n ~~  2 — 2 

\A“  - rn  COS- 


ocos"^o 


r 

= rb  sin  0b  - rQ  sin  fiQ  + f - XA 1- : 

J rK  VAr2  + Br  + < 


Again,  using  the  tables  of  standard  forms,  the  result  is 


, \/X-v/x^ 

Pr=rbsin<3b-r0sm(30+  A + 


B ?n  fv^V+Arh^' 
2 v/a^  n 1 >/**  + Ar  + B/2  . 


At  the  ray-path  peak,  where  r = r{,  equation  13  becomes 


p;  = rb  sin  0b  - rQ  sin 


The  total  group  distance  P'  = 2P^ . 


B n f\Ax^  + Arb  + B/2]  y/x^ 
- — Cn  4 - • - — r — 

ffiA  1 -\/(B/2)2-A  C . 


PHASE-PATH  DISTANCE 


The  phase  path  is  given  by  the  integral 


2 2 

Pr  = /Mds  = rbsin0b-rosin0o+  f 

JT  r vAr“  + Br  + 


where  r2p2  is  obtained  from  equation  4.  The  evaluation  of  the  integral  results  in  the 
equation 
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Pr  = rb  sin  0b  - rQ  sin  pQ  +y/x~y/x^ 

+ JL  Cn  2Ar  + B + 2n/AX  + ^_[m  gn  rb(2C  + Br  + 2y/CX) 
VA  2Arb  + B + 2v/AX^  2 y/C  r(2C  + Brb  + 2v/CX^) ' 


At  the  ray-path  peak,  r = rt,  equation  16  simplifies  to 
pt  = rb  sin  - r0  sin  _v^b 

- \/B“  - 4AC  , Br 


(16) 


+ — Cn 


m 


Cn 


rK  \/b2- 


4AC 


(17) 


2y/A  2Arb  + B + 2v/AX^  2v/C  2C  + Brb  + 2V/CX^ 

The  total  phase  path  is  P = 2Pt. 

MULTILAYER  MODEL 

During  the  night  hours,  the  ionosphere  is  usually  simple  enough  that  a single-layer 
model  is  adequate.  However,  during  the  daytime  there  are  two  or  three  layers  (E,  FI , and  F2). 
These  can  be  represented  by  separate  or  connected  parabolic  layer  segments.  During  the 
winter,  the  FI  layer  is  often  better  represented  by  a linear  layer  segment  (reference  2).  This 
can  be  accomplished  by  using  the  quasilinear  segment  (reference  3). 

The  equations  6,  13,  and  16  can  be  written  as  a sum  of  integrals  of  ray  properties  in 
each  layer,  where  the  upper  and  lower  integration  limits  are  the  altitudes  (radius  from  earth 
center)  of  the  layer  intersections.  The  equations  take  one  of  two  forms,  depending  on  whether 
( 1 ) the  ray  penetrates  all  the  layers  or  ( 2)  the  ray  is  reflected  in  one  of  the  layers.  Models 
having  layers  separated  by  free-space  regions  will  add  a slight  complication  which  will  be 
considered  later. 

Let  the  altitude  (radius)  of  the  ktB  layer  be  represented  by  a parabolic  layer  having 
coefficients  given  by  4,  so  that 


X=  Akr~  + Bkr  + Ck  , 


rk<r<rk+l 


To  ensure  that  the  electron  density  profile  be  continuous,  we  require  that  X have  the  same  value 
using  Ak,  Bk,  Ck  at  r = rk+|  as  using  Ak+j , Bk+j . Ck+j . Usually,  we  are  given  a model 
defined  by  equation  1 with  Nm,rb,and  ym  assigned.  The  boundaries  rk  are  roughly  known 
and  should  be  solved  for  using 


( Ak  - Ak-1  )rk  + ( Bk  _ Bk-1  )rk  + (Ck  “ ck-l  > “ 0 


(18) 


3.  Weast,  RC,  and  SM  Selby,  CRC  Handbook  of  Tables  for  Mathematics,  Chemical  Rubber  Co,  Cleveland, 
Ohio,  1970. 
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There  are  solutions  to  equation  18,  so  the  solution  consistent  with  the  model  is  to  be  deter- 
mined and  used.  It  is  possible  that  no  solution  is  consistent  with  the  model.  For  example, 
if  the  E and  FI  layers  are  modeled  by  some  world  map  function,  it  is  possible  that  under 
some  conditions  the  solutions  of  equation  18  will  not  be  in  the  ionosphere  (N  negative). 
These  conditions  must  be  checked  for  in  practical  applications.  We  now  consider  two  sets 
of  ray-path  equations.  First,  rays  which  completely  penetrate  the  ionosphere  and  second, 
rays  which  reflect  from  one  of  the  layers.  In  case  1 , we  consider  the  path  of  a ray  which 
penetrates  all  n layers  of  the  model  ionosphere. 


2 „ _J_  0 f rk+l(Ck+VCkXk  + Bkrk/2>1 

r-'o  b-  0>  r0‘°s  0 £ ^ n|rk,Ck+VEiA^tBk,k+|/2)J 


. (19) 


Pr  = rbsin0b~rOsin0O+  £ 

k 


t jVk 


+ \/Ai,  Xi.4.1  + Bu/2  J 


2^  Akrk+1  +v/AkXk+l  +V2- 


n 

Pr  = rbsin^b"TOsin%+  Z ■ ^Xk+1 

k=l  . 


Bk  Akrk+1  +VAkXk+l  + Bk/: 

Ak  rk+v/Akxk  + Bk/2 


(20) 


+ gn  rk<Ck+^Ck  Xk-H  +Bkrk+l/2>1 

2V^  rk+l(ck+v^ckxk  + Bkrk/2)  J 


(21) 


where  Xj  = Aj  r^  + Bj  rj  + Cj  and  r_  is  rm  for  the  kth  layer. 


In  case  2,  the  ray  reflects  in  the  layer  at  an  altitude  r^  where  X = 0.  Not  only  does  this 
simplify  the  equation  for  the  M**1  layer,  but  it  is  necessary  that  X = 0 exactly.  Numerical 
evaluation  of  equations  19  through  21 , using  equation  7 for  rj^+i , will  be  inaccurate  and 
result  in  square  roots  of  negative  (but  small)  numbers. 
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2 1 rk+l^Ck  +v/ck  xk  + Bk  rk^ 

Dr  = - Pn)  + rn  cos  pn  52  ~Cn — 

k=  i V^T  rR(Ck  + ^Ck  xk+l  + Bk  rk+ 1 l2) 


2 I 0 2(CM  +v'CMXM)  + BMrM 

+ r0  cos  Po~pzen 
VCM 


rM  K ~ 4AmC 


M'-  M 


r,  . , , V |v^7i 

Pr  = TbsmPb-r0smp0+  1,  1 JT 

k=l  ^ • 


Bk 

+ — Cn 


:rk+VAkXk  + Bk/2  1 
c+1  + \AkXk+1  + Bk/2  J 


2v^  Akrk+1  VAkXk+l  +Bk 


V^M  BM  „ 2<AMrM+v/AMXM)  + BM 
+ — — 8n  


Am  n/?"  -Vbm-4am  CM 


M 


M-l  • 

Pr  = rb  sin  0b  - rQ  sin  p0  + 52  ' V\+\ 

k=l  . 


Bk  „ Ak  rk+l  +VAkxk+l  + Bk^2 

+ — — Cn  


2V/Ak  Ak  rk +v/Akxk  + Bk/2 

en  rk^Ck  +v^kXk+l  * Bk  lW2| 
2^  " rk+|(Ck+v^  + Bkrk/2)  J 


-v^7+ 


B 


M 


-J3T- 


6n 


'M 


4AmC 


NTM 


M 2v/Aj^  2<AMrM  + v/AmXm)  + BM 
rM  " 4AMCM 


BMrm 


M 


en 


2VC^  2(cM  + v/CMXM)  + BM  rM 


(22) 


(23) 


1 


(24) 
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DISCONNECTED  LAYERS 


Sometimes  the  layers  are  disconnected,  which  results  in  a “valley”  between  layers. 
This  is  usually  the  case  at  dawn  between  the  E and  F layers.  This  adds  a complication  re- 
quiring special  treatment  of  the  region  between  the  upper  and  lower  layer  groups.  The  ray 
will  travel  in  free  space  between  the  layers.  Let  the  free  space  region  be  r,  < r < rj+  j . 
The  sums  in  equations  19  through  24  will  have  the  i^1  layer  deleted  and  replaced  by  a free- 
space  term.  For  Dr  in  equations  19  and  22,  use 


r0^i+l  “0i> 

for  the  i**1  term  in  the  sum.  For  Pj.  and  Pr,  use 
ri+l  sin0i+1  -rjsin^ 

instead  of  the  i**1  terms  in  equations  20,  21,  23,  and  24. 


SAMPLE  RAY  PLOTS 

The  equations  for  both  single-layer  and  multilayer  models  have  been  programmed  by 
means  of  minicomputers.  A sample  ray  plot  and  the  BASIC  program  which  produced  it  are 
displayed  in  figures  2 and  3.  Figure  4 shows  an  example  of  a three-layer  model  and  the  ray- 
trace  fan.  This  plot  was  produced  on  a Calcomp  plotter,  using  an  SEL-8 10  minicomputer. 


VALLEY  LAYERS 


The  region  between  layers  sometimes  has  a smoothly  varying  electron  density  with  a 
nonzero  minimum  in  N,  called  a valley  layer.  It  can  be  modeled  with  an  inverted  parabola  by 
d2N/dr2  > 0.  The  minimum  in  N may  be  zero  for  the  special  base  layer,  used  to  eliminate 
effects  of  the  discontinuity  in  N.  Since  the  signs  of  some  of  the  terms  have  been  changed, 
the  equations  for  P,  P\  and  D have  solutions  involving  the  arcsine  function.  Before  listing 
these  solutions,  we  rederive  the  constants  A,  B,  and  C for  two  cases. 


CASE  1,  BASE  LAYER,  Nm  = 0 


Since  Nm  is  unavailable  for  scaling  Ne,  we  use  Nb  defined  by  Ne  at  rb  where  rj,  = 
rm  " Ym  • *b 's  the  plasma  frequency  corresponding  to  Nb. 


Ne  = Nb 


ym  rj 
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Figure  3.  BASIC  program  used  to  produce  figure  2 on  a Tektronix  4051  microcomputer. 
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Figure  3.  (Continued). 
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and 


M2=l- 


f_b  rjn-r  r' 

f y 


m 


V ‘ 

r 


B “ 2 rm  [ ^W^m] 


C = -[fbrbrm/fym]2-rocos<30 


(25) 


CASE  2,  VALLEY  LAYER,  Nm  > 0 

The  constants  B and  C are  the  same  as  in  equation  25. 

21 


Ne  = Nm 


1 + 


r - rm  rb 


'm 


U2=  l-(fr/f)2- 


fb  rm  - r rb 


f Vm  r 


AM-cy^-CVb/fy,,,]2. 


(26) 


We  can  use  Snell’s  law  to  replace  r 2 cos^  p by  t2h2  cos2  P for  studies  with  P known  at  a 
certain  r. 

Only  some  of  the  terms  in  equations  19,  20,  and  21  are  changed  in  the  inverted  QP 
layer  solution.  Let  the  k*h  layer  (rk  < r < rj^  + j ) be  an  inverted  layer.  Then  the  following 
terms  in  each  sum  are  replaced.  In  equation  19,  use 


Bkrk  + 1 + 2(-k 


i 


- arcsin 


in  Bkrk  + 2Ck  j 
rk  ^i-4AkckJ 


(19i) 


In  equation  20,  replace  the  term  containing  Cn  with 


l 


2Akrk+l+Bk  2Ak 

--.arcsin  — =====  - arcsin  ~j= 

2v^Ajl  - yBJ-4AkCk  -v/B l 


rk  + Bk.  1 
-4AkCkJ 


(20i) 
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And  replace  the  last  term  in  equation  2 1 with 


' 


; 


I 


i 


V-^k 


®krk+l+^k  ®krk  + ^k 

arcsin  arcsin 


*k+l^-4AkCk 


rkyB^-4AkCk 


Figure  5 shows  an  example  of  a multilayer  profile  with  inverted  segments. 


Figure  5.  Multilayer  model  (not  to  scale)  indicating  the 
three  different  types  of  quasiparabolic  forms:  normal  QP 
layer;  valley  layer;  and  base  layer  used  for  smooth  transi- 
tion to  free  space,  with  no  discontinuity  in  the  slope. 


(21i) 


If  a transmitter  were  located  in  the  valley  layer  (r3  < r < r4>,  it  would  be  possible  to 
launch  rays  which  are  trapped  between  the  upper  and  lower  sides.  The  hop  period,  Dj,,  is 
obtained  by  substituting  equation  7 for  r^+i  and  r^  in  equations  19i  and  19.  Equation  7 can 
be  written  as 
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■ -IBBWP 


2Art  + B = B + 2C/rt  = ± \/b2  - 4AC  . 

Choosing  + for  r^+j  and  - for  in  equation  19i,  we  have 


Dh  = 


2irrQc°sP0 
T- T~ = 2mn 

y/~^k 


2 2 a / 2 2 a ~‘m 

r0COS  V r0cos  0o  - — 


Br, 


1/2 


Using  Snell’s  law,  r2  cos20  = r2  n 2 cos20,  we  obtain 

Dh  = 2tttq  |l  - s 111 2~l 

u ^ 2Ar2  + B(2r  + rm)  cos^/jJ 


where  (3  is  the  angle  of  the  ray  at  altitude  r.  The  trapped  ray  travels  in  a wavelike  path  with 
period  Dj,  and  peak-to-base  altitude  range  5r  given  by 


6r  = 


» 

1 

— « 
A 


Bz  - 2ABrm  + 


4A2  r2  + 2AB(2r  + rm) 


0 -1 


c«4l/2 


If  the  altitude  range  is  small  and  restricted  to  the  upper  portion  of  the  valley,  we  have 
a “whispering  gallery”  mode.  For  0 = 0,  we  get  r = 0 when  r = Tg,  the  whispering  gallery 
radius: 


r --Bk/2Ak  . 


SUMMARY 

The  quasiparabolic  layer  ray-tracing  equations  first  presented  in  reference  1 for  a 
single  layer  have  been  extended  to  a multilayer  model  including  valley  layers.  The  equa- 
tions for  points  on  the  ray  trajectory  are  used  to  display  ray  paths  reflecting  from  and 
traversing  through  a multilayer  ionosphere.  The  parameters  of  the  whispering  gallery  rays 
in  the  valley  layer  are  presented. 

It  is  hoped  that  these  equations  will  find  application  in  multifrequency  communi- 
cation networks.  The  equations  presented  are  compact  enough  so  that  they  can  be  solved 
by  means  of  mini-  and  microcomputers  in  an  interactive  manner. 
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